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Scientists studying the atmosphere typically rely on mathematical and computer mod- 
els to try to make sound predictions about weather and climate change. By themselves, 
these models are not enough to allow for very accurate predictions since not all natural 
processes are known and accounted for in the models. Data assimilation is the vehicle 
used by scientists to bring model predictions close to reality. Atmospheric data assimila- 
tion consists of a series of mathematical steps that combine model predictions with actual 
observations of the atmosphere to produce an estimate of the state of the atmosphere at 
any given time. The estimates are commonly referred to as analyses. When all goes well, 
the analysis is a better estimate of the state of the atmosphere than the estimate provided 
by either the model or the observations alone. Usually, only observations before and at the 
time of the analysis are used to calculate this “filter” estimate. The name filter comes es- 
sentially from the fact that, in a manner of speaking, this assimilation procedure combines 
the best of two worlds by filtering out their errors: the “observation-only world” and the 
“model-only world” . More sophisticated assimilation procedures known as smoothers are 
capable of combining filter estimates with observations within a certain time interval to 
produce refined estimates of the state of the atmosphere, within the desired time interval. 

There are different smoother types. In the present work the so-called fixed-lag Kalman 
smoother is used as a framework to construct a retrospective assimilation system for the 
NASA/Goddard Earth Observing System (GEOS) Data Assimilation System (DAS). In 
this type of smoother formulation, when observations up to 6 hours ahead of a regular filter 
estimate are used to calculate the (refined) retrospective estimate we say we are calculating 
the lag-1 retrospective analysis; when observations up to 12 hours ahead of a regular 
filter estimate are used to calculate another (even more refined) retrospective estimate 
we say we are calculating the lag-2 retrospective analysis; and so on. The results of our 
experiments with GEOS DAS indicate that the lag-1 retrospective assimilation procedure 
does indeed provide an overall improvement over the regular assimilation procedure. One 
particular significant result, obtained by studying the skill of 5-day forecasts, indicates 
that lag-1 retrospective analyses seem to consist of better initial conditions than those 
normally provided by the filter analyses. Even though our results are obtained for a 
slightly simplified version of GEOS DAS, they are quite promising and work is already in 
progress to expand this research, including study of the impact of lags higher than one. 
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Abstract 


The fixed-lag Kalman smoother (FLKS) has been proposed as a framework to construct data 
assimilation procedures capable of producing high-quality climate research datasets. Fixed-lag 
Kalman smoother-based systems, referred to as retrospective data assimilation systems, are an 
extension to three-dimensional filtering procedures with the added capability of incorporating 
observations not only in the past and present time of the estimate, but also at future times. A 
variety of simplifications are necessary to render retrospective assimilation procedures practical. 

In this article, we present an FLKS-based retrospective data assimilation system implemen- 
tation for the Goddard Earth Observing System (GEOS) Data Assimilation System (DAS) . The 
practicality of this implementation comes from the practicality of its underlying (filter) analysis 
system, i.e., the physical-space statistical analysis system (PSAS). The behavior of two schemes 
is studied here. The first retrospective analysis (RA) scheme is designed, simply to update the 
regular PSAS analyses with observations available at times ahead of the regular analysis times. 
Although our GEOS DAS implementation is general, results are only presented for when obser- 
vations 6-hours ahead of the analysis time are used to update the PSAS analyses and thereby to 
calculate the so-called lag-1 retrospective analyses. Consistency tests for this RA scheme show 
that the lag-1 retrospective analyses indeed have better 6-hour predictive skills than the pre- 
dictions from the regular analyses. This motivates the introduction of the second retrospective 
analysis scheme which, at each analysis time, uses the 6-hour retrospective analysis to replace 
the first-guess normally used in the PSAS analysis, and therefore allows the calculation of a 
revised (filter) PSAS analysis. Since in this scheme the lag-1 retrospective analyses influence 
the filter results, this procedure is referred to as the retrospective-based iterated analysis (RIA) 
scheme. Results from the RIA scheme indicate its potential for improving the overall quality of 
the assimilation. 
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1 Introduction 


The concept of retrospective data assimilation, as invoked in the present article, was introduced 
by Cohn et al. (1994; CST94 hereafter) to refer to the calculation of the analyses from observa- 
tions after the analysis time, as well as before and at the analysis time as is done in numerical 
weather prediction. Retrospective data assimilation is possible when analyses are not required 
in real time, such as in the production of reanalysis data sets for climate research. 

In estimation theory, estimates of the state of a system produced from observations on both 
sides of the analysis time are known as smoother estimates. In sequential data assimilation a 
natural smoothing technique to employ is that of fixed-point smoothing. In this case, the usual 
filter estimate obtained at a fixed time using observations before and at the analysis time is 
sequentially updated as future observations became available. Future observations can be used 
for as long as experimentation shows their impact to be useful. The idea of estimating the 
state of a system at a fixed time over and over again as more observations become available can 
be taken a step further by seeking fixed-point estimates at a series of consecutive fixed times. 
This is what is accomplished by fixed-lag smoothing. Specifically, for linear systems under 
the typical assumption of unbiased Gaussian-distributed errors the fixed-lag Kalman smoother 
(FLKS) provides the best unbiased estimate of the state of the system at a sequence of given 
times using observations in the past, present, and at a time lag-^ ahead of the time of each 
estimate. 

The FLKS is composed of two major components: the Kalman filter (KF) portion and the 
fixed-lag smoother portion. The FLKS is fully dependent on the KF as it is formulated on the 
basis of the observation-minus-forecast residuals resulting from the KF. In general, when the 
filter is not the KF, but rather some suboptimal implementation of it, we can still think of sub- 
optimal implementations of FLKS-based retrospective data assimilation schemes as consisting 
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of a filter portion and a smoother (or retrospective) portion. Todling et al. (1998) used this 
explicit separation between filtering and smoothing portions to study the behavior of a vari- 
ety of combinations of filter and smoother approximations to the linear FLKS. One particular 
approximation studied there, namely the adaptive CCF-based retrospective data assimilation 
scheme, was seen as having the potential for being implemented in practice. It replaces the 
filter portion of the FLKS by a constant forecast error covariance filter much like operational 
three-dimensional variational analysis systems do. The spectral statistical-interpolation analysis 
system of Parrish and Derber (1992) is an example of such a system; the U.S. Navy analysis 
system of Daley and Barker (2001) is another; the European Center for Medium-Range Weather 
Forecasts (ECMWF) system of Courtier et al. (1998) is yet another; and so is the physical-space 
statistical analysis system (PSAS) of Cohn et al. (1998), which is also central to the work in 
the present article. 

To take forward the idea of developing a practical retrospective data assimilation system, 
the linear FLKS formulation of CST94 has to be extended to handle nonlinear dynamics. Since 
the retrospective portion of the algorithm relies completely on the filter, designing nonlinear 
filters immediately results in designing nonlinear smoothers. Todling and Cohn (1996; TC96 
hereafter) derived a nonlinear FLKS algorithm based on the traditional extended Kalman filter 
(EKF). Similar derivations can be found elsewhere (e.g., Biswas and Mahalanabis 1973; Verlaan 
1998). The way smoothers use future observations to calculate updates to state estimates, is 
by propagating information back in time using the adjoint dynamical model. For nonlinear 
dynamics the adjoint of the tangent linear dynamics must be provided in principle. Four- 
dimensional variational (4D-var) procedures such as that of Rabier et al. (2000) also require the 
adjoint of the tangent linear dynamics. The need for the adjoint model can be avoided if the 
retrospective assimilation strategy is based on ensemble techniques such as that of Evensen and 
van Leeuwen (2000). 
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In this article, we study the performance of a PSAS-based retrospective analysis (RA) system 
developed for the Goddard Earth Observing System (GEOS) Data Assimilation System (DAS). 
Since the forecast error covariance matrix of PSAS in GEOS DAS varies slowly in time we 
can identify the suboptimal RA procedure studied here with the CCF scheme of Todling et al. 
(1998). Our RA implementation in GEOS DAS is general and applicable to any number of time 
lags, but in the present article we concentrate on results for the 6-hour, i.e. lag-1, retrospective 
analysis. Motivated by some of the results obtained with this version, and by the ideas of 
constructing so-called iterated filters and smoothers common in the engineering literature, we 
also study here the performance of a retrospective-based iterated analysis (RIA) scheme. In the 
RIA, the lag-1 retrospective analysis at a given time tk-i is used to produce a new first-guess at 
time tk that is used to revise the filter (PSAS) analysis at the same time £/-. In the RIA the final 
analysis is the second (iterated) analysis calculated using the first-guess generated from the lag-1 
retrospective analysis. This is a considerably different use of the “static” retrospective analyses 
proposed by CST94. Though a formal argument for the RIA procedure is not presented here, 
the procedure is found to improve the overall quality of the analyses. This lag-1 RIA scheme 
makes the retrospective procedure resemble a 4D-var cycle (e.g., Courtier et al. 1994; Rabier et 
al. 2000, Li and Navon 2001). 

Indeed, the original FLKS-based retrospective analysis formulation of CST94, and the RIA 
here, can be viewed as alternative approaches to 4D-var. The FLKS framework is a natural 
four-dimensional extension to three-dimensional procedures formulated sequentially rather than 
variationally. Four-dimensional variational procedures are an extension of 3D-var that take 
into account observations within a time interval. Menard and Daley (1996) have shown the 
equivalence of 4D-var and fixed-interval smoothing. Similarly, for linear dynamics, the FLKS is 
algebraically equivalent to 4D-var and can be derived from the 4D-var cost function by solving a 
two-point boundary value problem (Zhu et al. 1999). The main distinction between 4D-var and 
the FLKS is in their computational approaches. The former involves an iterative optimization 
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procedure to arrive at the solution, whereas the latter deals directly with the analytical solution 
of the problem. One practical consequence of this distinction relates to how these procedures 
account for model error. As pointed out by Todling et al. (1998), FLKS-based assimilation 
schemes directly inherit any model error covariance parameterization embedded in the filter 
portion. Various techniques to account for model error in 4D-var can be formulated by using the 
dynamical model as a weak constraint on the optimization problem (e.g., Derber 1989; Bennett 
et al. 1996; and Zupanski 1997). However, until a more complete understanding of model error 
is acquired, and the corresponding model error covariance parameterizations can be relied upon, 
this distinction between 4D-var and FLKS-based assimilation is rather moot. Another important 
point to make relates to what now seems to be recognized (Fisher and Andersson 2001) as one 
of the main advantages of 4D-var over 3D-var-like procedures, namely, that the former uses the 
observations nearly at their proper times [as the case of the ECMWF 4D-var implementation of 
Rabier et al. (2000)], whereas in the latter it is more common to bundle the observations into 
6-hour batches. This can be resolved, particularly in sequential 3D-var assimilation procedures, 
by using a rapid update cycle strategy. Though this is not explored in the present article, since 
in GEOS DAS the observations are bundled into 6-hour batches, we should point out that there 
is no intrinsic difficulty in building an FLKS-based retrospective analysis system under a rapid 
update cycle filtering strategy. 

In the sequel we briefly review, in section 2, the theoretical framework behind retrospective 
analysis. The presentation is based on the EKF and the corresponding nonlinear extension of 
the FLKS. >In section 3, we describe the framework of our practical implementation directed 
toward adding a retrospective component to GEOS DAS; here, both the RA and RIA schemes 
are presented. In section 4, results of a preliminary evaluation of these retrospective schemes 
are presented and discussed. Conclusions are drawn in section 5. 
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2 Theoretical framework: the fixed-lag Kalman smoother 

In this section we briefly recapitulate the formulations of the fixed-lag Kalman smoother of 
CST94 and TC96. Following Todling et al. (1998) we separate the FLKS into a filter portion and 
a retrospective portion. The filter portion is based on the linear Kalman filter, or more generally 
on any nonlinear extension of the KF; the retrospective portion is based on the linear fixed-lag 
Kalman smoother, or any equivalent nonlinear extension compatible with the underlying filter. 
As in TC96, the discussion below is based on the EKF. 

(a) The filter portion 

Using the notation of CST94, the filter portion of the FLKS formulation of TC96 can be 
summarized by the usual EKF equations 


k\k-l 

= » 

(la) 

w £i* 

= W £|/c-l + K k\k V k , 

(lb) 

^-k\k 

= r Vi H * r : 1 . 

(1c) 

i 

k\k— 1 

= A-k,k-l'Pk-l\k-l-^-k,k-l + > 

(Id) 

pa 

r A:|A; 


(le) 


The first two expressions refer to the state estimate evolution, which depends on the last three 
expressions essentially related to error covariance evolution and update. At time tk, the forecast 
72-vector w^_ 1 evolves through the nonlinear dynamical operator Ak,k - 1 from the analysis 72 - 
vector w^_ 1 | A; _ 1 , according to (la). The dynamical operator Ak,k-i stands for, say, a general 
circulation model, and possibly any transformations. necessary to convert the model prognostic 
variables into the filter state vector, and vice-versa. 

The main difference in the EKF equations written above and the way they more commonly 
appear in the atmospheric data assimilation literature (e.g., Miller et al. 1994) is in the time 
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subscripts. Here, the subscripts follow standard engineering notation developed in estimation 
theory and which is mostly suitable to the development of smoothers. This subscript notation 
is also particularly helpful in reminding us that for linear systems perturbed by Gaussian- 
distributed noise the forecast w^_ 1 and analysis state vectors are actually conditional 
means of the true state n- vector wjj., that is, 

w ^-i = £{ w fci w fc-n--- >o> ( 2a ) 

w t\k = > w ?}> ( 2b ) 

at time t &. The conditioning, represented by the vertical bar in the expectation operator £{• |«}, 
is on the time series of observations w£. The forecast at time tk is the expected value of the true 
state conditioned on all observations prior to time tk] the analysis at time tk is the expected 
value of the true state conditioned on all observations up to and including those at time tk . 

The EKF, like the KF, depends on the residual p*;-vector Vk in (lb) formed by the difference 
between the vector of observations w£ and the model-predicted “observations” 
at time £&, that is, 

= w£ - ?4(W f k]k _ x ) . (3) 

The nonlinear observation operator Hk stands for the transformations involved in converting 
filter state vector quantities into observables. Optimality of the filter depends on the n X 
Pk weighting matrix given to this observation-minus-forecast (OMF) residual vector 

through (lb). Although the expression for the weighting matrix K k \k in the EKF is similar in 
form to its linear KF equivalent, contrary to the linear case, in (lc) is now state-dependent 
since the pk X n Jacobian matrix of the observation operator Hk is linearized around the 
- forecast state vector w State dependence of the EKF weighting matrix also comes 
from its dependence on the OMF residuals covariance matrix T^, given by 

r fc = H fc P{| fc _ 1 H|’ + R fc , (4) 
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for uncorrelated observation and forecast errors. Here, is the pk X pk observation error 
covariance matrix and is the state-dependent nxn forecast error covariance matrix. The 

dependence of the forecast error covariance matrix in (Id) on the state comes from the nxn 
Jacobian matrix Ak,k-i of the dynamics operator Ak,k - 1 which is linearized about the model 
trajectory initialized from the analysis vector w^_ 1 | A ._ 1 . The forecast error covariance matrix 
also depends on the model error covariance matrix Q/~, which is normally assumed to be known. 
Since the forecast error covariance matrix P{|^_ 1 evolves from the nxn analysis error covariance 
matrix P^_ 1 | A: _ 1 it depends further on the accuracy of the previous estimate calculated by the 
filter, i.e., through (le) applied at time tk- 1 - 

In the linear case, the dynamics and observation operators reduce to Ak,k - 1 = ^k,k - 1 and 
'Hk = H&, respectively, and (1) reduces to the linear KF for known model and observation error 
statistics. Moreover, as pointed out in TC96, in the linear Gaussian-distributed noise case, the 
forecast and analysis error covariance matrices are the conditional mean error covariances. It 
is when the observation errors are Gaussian and white in time, that the time series of residual 
vectors v& can be identified with the innovation sequence (see for example, Anderson and Moore 
1979, section 5.3) 

(b ) The retrospective portion 

In the FLKS, the retrospective portion uses the OMF residual vector at time tk to 
calculate corrections to filter analyses and retrospective analyses at previous times tk-i using 
an update equation similar to the state update expression (lb) of the filter portion. The lag-^ 
FLKS retrospective analyses based on observations newly available at time tk are calculated by 

™l-i\k = w k-t\k-i + , (5) 

for £=1,2,..., min(&, L), and a maximum desired lag l = L. They are analyses for times tk-i- 
Each retrospective analysis for fixed time tk-i is also an “incremental” correction to an estimate 
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of the state calculated previously. For example, when k = k and t = 1, the lag-1 retrospective 
analysis is a correction to the most recently available state estimate at time t K - 1 , i.e., the 

filter analysis w“_ 1 | K _ 1 , based on the observations newly available at time t K \ when k = k + 1 
and l — 2, the lag-2 retrospective analysis w“__ 1 | /t+1 is a correction to the most recently available 
state estimate at time t K - 1 which is now the lag-1 retrospective analysis w^_ 1 | k ; and so on up to 
the desired lag i — L when the estimate at time t n -i is given by the lag-L retrospective analysis 

w a 

K — — 1 

This example to illustrate the mechanism for correcting consecutive state estimates at a given 
time with successive smoother calculations makes the FLKS algorithm resemble very much the 
fixed-point smoother. This is simply because in this example we chose to fix the time at which 
estimates are being sought, that is, time i. The resemblance between the fixed-lag and 
fixed-point smoothers is no coincidence. The FLKS of CST94 and TC96 can be derived from a 
fixed-point smoother formulation using, for example, the approach of state augmentation (e.g., 
Biswas and Mahalanabis 1973). Out point here is simply that the incremental corrections to 
the state estimates at a fixed time are calculated on the basis of the OMF residual vectors 
v«, v« +1 , and so on up to v*+£,_i. That is, each lag of the algorithm introduces corrections to 
the state estimate by using observations at times further and further ahead of the retrospective 
analysis time, up to the maximum desired lag L. 

Because the retrospective analyses are based on the same OMF residual vectors used in 
the filter portion of the algorithm, the retrospective n X pk weighting matrix K k-i\k depends 
on the OMF residual covariance matrix in (4). Furthermore, K^i*. also depends on the 
n X pk matrix , the transpose of the Jacobian of the observation operator, and on the n X n 
forecast-analysis cross-covariance matrix through the EKF-based expression 

K i _ (lt = (P*_ (|t _ I ) T Hi'rj 1 , (6) 

as can be found in TC96. The forecast-analysis cross-covariance | i _ 1 evolves from previ- 
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ously calculated analysis error covariances and analysis-analysis error cross-covariances through 
the Jacobian Ak,k-i of the dynamics operator. Its evolution equation and the update equations 
for the retrospective analysis error cross-covariances are 

P k-i\k — P k-t\k-i - '^k-e\k'KkP{ tk _£\k-i > ( 7a ) 

= (I - , (7b) 

^l,k-t\k-i ~ A-W-i'P k-i,k-e\k-i » 

and the details of their derivation can also be found in TC96. 


That retrospective analyses are built on the basis of future observations can be simply 
understood by recalling the meaning of the time subscript notation used here. In the linear 
Gaussian-distributed noise case the time subscript notation signifies that the retrospective anal- 
ysis estimates are indeed estimates of the conditional means. In this case, the retrospective 
analysis at time tk-i is 


™l-i\k = £Wk-i\™h w L i, • • • . O > 


( 8 ) 


where now, in contrast to the filter estimates (2) , the expectation is conditioned on all obser- 
vations before, at and after time U P to time ffc. As mentioned previously, in the linear 
optimal case, when the underlying filter is the KF and the sequence of OMF residual vectors is 
actually the innovation sequence, the retrospective portion just described reduces to the optimal 
FLKS. Independently of nonlinearities, in general, if the filter is suboptimal the corresponding 
retrospective analyses are suboptimal as well. This is simply because both the filter and the 
smoother are based on the same sequence of OMF residual vectors v*,. Unfortunately, in the 
suboptimal case, there is no guarantee that consecutive retrospective lagged estimates will repre- 
sent improvements over estimates with smaller lag(s) or even over the filter results (see Todling 
et al. 1998 for illustration) . 


As pointed out by Todling et al. (1998), one interesting feature of the FLKS that arises 
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directly ton its being formulated on the basis of an undoing «Her in that it incorporates model 
error covariances naturally and automatically (see also appendix A, here). In fact, equations 
( 5 )-(7) do not depend explicitly on the model error covariance. A variety of techniques exist 
t o incorporate mode, error in 4D-var (e.g„ Derber 4989; Bennett et a, 4996; ^ Znpanshi 
W97). Since 4D-var is algebraically equivalent to fixed-interval smoothing (see Menard an 
Daley 1996; and Zhn et al. 1999) and for all practical purposes we can always choose a lag L 
in fixed-lag smoothing that accomplishes the same beneht as ^-interval smoothing (Moore 
W73 ), FLKS- based assimilation procedures present a potential aiternative to 4D-var. Since we 
currently lack the necessary knowledge to parameterize model error covariances this advantage 
of the FLKS over 4D-var is not very significant, but it may prove to be relevant in the future. 

3 Practical framework: GEOS DAS considerations 

The aigorithm described in the previous section serves mainly as a guide to help design suitably 
feasible data assimilation procedures. It is well known that the computational cost of evo vmg 
full covariances is excessive for filtering, let alone for smoothing as in (7), and likeiy not yustifia e 
because of our reiative lack of knowledge of the required input model and observation error 
statistics. This has motivated the study of a number of simplihcations to both hltermg (e.g 

Cohn and Todling !996, and references therein) and smoothing (e.g., Todling et al. 4998, and 

. , Q v this section we describe the details of our implementation o 

references therein) procedures. In. this sectio , 

j t f u nFOq DAS Before describing the retrospective 
the FLKS-based retrospective procedure for the GEOS DAS. tfetor 

• current GEOS DAS that approximates, 

analysis portion of the implementation we summarize the 

in principle, the filter portion of the algorithm. 

( a ) The GEOS analysis and data assimilation system 

The DAO operational GEOS data assimilation system consists of three major comp 
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an atmospheric general circulation model (GCM) ; the physical-space statistical analysis system 
(PSAS); and the incremental analysis update (IAU) procedure. At the so-called analysis times, 
the GCM provides a first-guess field to PSAS so it can process OMF residuals and generate the 
analysis state. The physical-space statistical analysis system is an implementation of the EKF 
equations (lb)-(lc), obtaining the analysis state as a correction to the model first-guess. The 
error covariance evolution expressions (Id) and (le) are neglected and therefore PSAS functions 
as a suboptimal filter, as in the case for other operational 3D-var systems. Each PSAS analysis 
is used in the IAU procedure of Bloom et al. (1996) to construct a tendency term that is used to 
force the GCM during a 6-hour period around the analysis time. The GCM trajectory obtained 
during the IAU integration is known as the assimilated trajectory. 

In GEOS DAS the state-space of the GCM is different than the state-space of the analysis 
system and it is convenient to define a specific nomenclature for the purposes of the present 
article. In what follows, we refer to background as the state-vector provided by the GCM and 
to forecast or first-guess as the background field transformed to the analysis space. The model 
and analysis spaces are different because their state variables and grids are different. The GCM 
state variables are surface pressure, potential temperature, specific humidity and the zonal and 
meridional components of the wind, where all variables are defined on the Arakawa C-grid and 
on a vertical sigma coordinate system. On the other hand, the analysis state vector is composed 
of sea level pressure, the zonal and meridional components of the sea level wind, the zonal and 
meridional components of the upper-air wind, mixing ratio, and geopotential heights, where all 
variables are defined on the Arakawa A-grid and in pressure coordinates (see DAO 1996, for 
details) . 

We designate an m-dimensional sigma-coordinate GCM state vector by y (o') and an n- 
dimensional pressure-coordinate analysis state vector by w(p), to emphasize explicitly the ver- 
tical coordinate system these states are defined on. For our purposes, we can represent a GCM 
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integration as 

- MyO)] + a 8 yl\k{ a ) •, ( 9 ) 

Here, M is the nonlinear GCM operator and the second term on the right-hand side corresponds 
to the constant IAU forcing term applied to the GCM during the IAU integration period. The 
parameter a controls when and how the model-space analysis increment aifecte the 

integrations. For 6-hours the IAU time interval [tk-1/21 Uc+1/2] we se ^ T%au = ^fc+1/2 ~~^k- 1/2 an< i 
a = 1 / T iau , and during the 3-hour GCM background integration time interval [t k +i/ 2 i t k+i] we 
set a — 0. At an analysis time tk, the GCM-provided background held y^-iO 7 ) is converted 
into the analysis hrst-guess through the operation 

w l\ k ~i ( p) = ( 10 ) 

where for convenience we use similar time subscript notation as that used in the previous section. 
The space conversion operator II is nonlinear since it represents not only simple interpolation 
from one grid to another but also variable transformations such as conversion from potential 
temperature to geopotential heights. This operator can be absorbed in the definition of the 
state vector and become transperent in the description of the filter and smoother equations. 
However, to make clear the connection between the mathematical description and the actual 
implementation of these procedures we opt to refer to n explicitly. 

The forecast vector w{| A ._ 1 (p) is used to construct the OMF residual p-vector in (3). 
Instead of calculating explicitly the weighting matrix (lc), PSAS splits the calculation of the 
last term in the analysis equation (lb) into two steps. The first step is to solve the linear system 
of equations 

r fc x fc = v fc , ( n ) 

for the variable X&, so that in a second step the analysis (p) can be calculated by 

w £|fc(p) = w fc|*-i(P) + P {|fc-i H fc x ^- ( 12 ) 
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To keep notation simple, we denote the PSAS forecast error covariance with the same symbol 
p / i US ed in the previous section. However, as mentioned above, PSAS does not use (Id) to 
calculate the forecast error covariance matrix. Instead, the forecast error covariance in PSAS is 
parameterized using simple dynamical constraints. Only its variance fields vary (slowly) in time, 
its correlations are constant in time. A consequence of such simplification is that the forecast 
w£, fc _i(p) and the analysis w^ fc (p) vectors in (12) are also distinct from those of the previous 
section, even though they are designated with the same symbols as in the previous section. 
Furthermore the forecast error covariance formulation of PSAS is for the analysis variables and, 
in particular, in pressure coordinates. Moreover, the observation operator Tik m PSAS is linear, 

that is, %k = Hfc. 

To proceed with the GEOS IAU assimilation, the analysis in (12) is converted back to the 
model space, through a conversion operator II + , 

y = n + [w£ |fe (p)] , ( 13 ) 

which is then used finally to construct the IAU 6y a k{k {a) increment to be used in (9), 

£y£ifcO) = yk\k(°) - y b k\k-i(?y • ^ 

The actual implementation of n+ is such that it renders minimal the difference between a field 
w(p) in the analysis space and the field resulting from transforming w (p) to the model space 
using H+ and subsequently transforming the resulting vector back to the analysis space using 

n. 

A schematic representation of the IAU assimilation procedure is shown m Fig. 1. In GEOS 
DAS observations are processed in 6-hour intervals, which in the IAU framework implies that 
the GCM is integrated for 6 hours starting 3 hours before the analysis time. Going from the 
left to right in the diagram, at an analysis time, say t = 6Z, observations and a 3-hour model 
first-guess (represented by the north-eastward pointing dashed arrow) are combined in PSAS 
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to calculate the filter analysis. This analysis is used to construct the IAU increment (14) and 
the model is integrated forward forced by the IAU tendency starting from t = 3Z up to t — 9Z. 
Beyond this time, the IAU forcing is set to zero and the model runs “free” for the next 3 hours. 
At the end of this free 3-hour integration the GCM provides the background to be used in the 
PSAS analysis of the 12Z observations, and the cycle is repeated. The assimilated trajectory is 
represented in the figure by the thick-solid eastward-pointing arrows. 

(b) The GEOS retrospective analysis 


We now have the challenge of converting the retrospective portion of the FLKS as presented 
in the previous section into a practical algorithm. We have seen above that when building a 
practical filtering procedure such as PSAS one of the main approximations is to avoid dealing 
directly with the error covariance equations (ld)-(le). Analogously, when building a practical 
implementation of the retrospective portion of the FLKS we want to calculate retrospective 
increments 


ti w k-i\k(p) — w k-t\k(p) ~ w k-i\k-i(p) — ^k-t\k^k » (!5) 


for lags £ = 1,2, ••• , mm(fc,L), without having to calculate the smoother cross-covariances 
implicit in the retrospective gains through (6) and (7). As it turns out, calculating these 

cross-covariances can be avoided since the retrospective gain matrices ~Kk-i\k can be written as 
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(16) 


Ij-k-i+l 


(see appendix A), with the consequence that the retrospective increments in (15) become 
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Hjfxfc, (17) 


\j=k-£+l 


where we used (11) to replace with X&. We see from this expression that the lag-£ ret- 

rospective increment is a linear combination of the columns of the forecast error covariance 
as original filter increment. The advantage of the expression above is that 



it refers only to quantities used by the filtering portion of the FLKS: the (filter) forecast error 
covariance matrix Pj_ 1 | j _ 2 ; the observation error covariance matrix Rj-i; the linear (or lin- 
earized) observation operator Hj_i and its transpose (adjoint); and the adjoint of the Jacobian 
Ajj_i of the dynamics operator. The smoother error cross-covariances Pk,k-t\k-\ an< l ^ k,k-i\k' 
and smoother error covariance P k -i\k do n °t appear in (17). 

At a given analysis time t k , the retrospective increments can be calculated through a succes- 
sion of operations similar to the two-step PSAS operations (11) and (12). Defining an n-vector 

z k\k as 

Zfclfc =Hjx fc , (18) 

corresponding to the PSAS conjugate gradient solution x fc converted from the observation space 
to the analysis space by HjT, the term in the square brackets of (17) can be calculated using the 
following algorithm: 


j = k 


while j > 1 and j > max(l, k — i + 1) 


z j-i| k = ^J,3- i z iP 

(19a) 

= H i-i p i- 11 i-2 z i- 1 |* 

(19b) 

v tj T r 

z j-l\k — z j- i|fc. 

(19c) 

5v?j-i\k(p) = Pj_i| J _2 z i-ip 

(19d) 


3=3 ~ 1 

endwhile 

for a maximum number of time lags l — L. In this algorithm the n-vector in (19a)is 

the result of the adjoint dynamics evolution of the auxiliary n-vector Zj\ k , for each backward 



integration j. This backward-propagated vector serves as the input to an equation (19b) 

similar to the first step (11) of the regular PSAS analysis, but now with a different right-hand 
side. The next step in the retrospective analysis loop is to update the n - vector with 

the analysis-space projection of ^ in (19c). Finally, the n - vector Zj-i\k i* 1 (19c) is used to 
calculate the retrospective analysis increment for each desired lag £ up to a maximum lag £ = L 
through application of the forecast error covariance operator in (19d). 


Notice that the entire retrospective analysis algorithm (18)-(19) works in the analysis space. 
In particular, the propagation operator Aj fc _ 1 = Aj A ,_ 1 (p) in (19a) is defined in pressure 
coordinates and it operates on geopotential heights, mixing ratio, zonal and meridional winds, 
etc, that is, the analysis variables. In fact, the linearized dynamical operator A^fc-i (p) is given 

by 


A~k,k-i{p) = n fc M fc) fc_i(a)n^_ 1 , 

where M fc ,fc_i(o-) is the m X m Jacobian matrix of the nonlinear operator M in (9), 


and II and II + are given by 


M(ct) = 

n = 

n+ = 


dM[ y] 


dy 


y=y{<?) 


dU[y] 


dy 
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( 20 ) 

( 21 ) 

(22a) 

(22b) 


W— w (p) 


and correspond to the n x m and m X n Jacobian matrices of II and II + , respectively, where 
we recall that m is the dimension of a model state vector and n is the dimension of an analysis 
state vector. 


A few remarks can be made at this point. 


• Currently in PSAS the analysis error covariance matrix is never referenced. Indeed, 
the current implementation of PSAS parameterizes the forecast error covariance matrix in 



such a simple manner that none of the terms on the right-hand side of (Id) are taken into 
account. However, when the expressions (7) for the smoother error cross-covariances are 
bypassed and the retrospective increments are calculated using the gains in (16) there are 
actually no approximations involved. The only consequence of not calculating the smoother 
error covariances is that we get no estimates for the accuracy of the retrospective analyses 
— which, in principle, can be extracted from P“_^. Expression (16) is exact for the linear 
FLKS and its nonlinear EKF-based extension. 

• We see from (17) that an FLKS-based retrospective scheme allows future observations to 
be used to correct previous filter and retrospective analyses impaired by the lack of obser- 
vations over a particular region earlier on in the assimilation. That is, when at time ffc-i, 
say, there are no observations over a certain region, the filter analysis at this time will 
essentially equal the first-guess over that region — aside from possible contributions by 
farther away regions through the forecast error correlations. If at time ffc, say, observations 
become available over the region in question, or information from observations at nearby 
downstream regions get propagated through the adjoint of the tangent linear dynamics 
A^k - 1 into t ^ ie region in question, this new information will be used to calculate a cor- 
rection to the filter analysis at time tk-i as the lag-1 retrospective analysis represented in 
(17). In these cases, it is the first term in the square bracket of (17) that mostly contributes 
to the correction to the filter analysis. 

• Notice that the linear system (19b) solved within the retrospective analysis algorithm 
involves exactly the same operators required to calculate the sensitivity of forecasts to 
observation changes, as measured by some pre-specified cost function, as in the approach 
of Baker and Daley [2000; compare with their eq. (2.7a)]. Furthermore, (19c) involves 
exactly the operator required to examine forecast sensitivity with respect to changes in 
the background. It has been pointed out elsewhere that some of the operations in 4D-var 
are closely related to operations required to study forecast sensitivity; the same is true of 
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the operations in FLKS-based retrospective analysis schemes. 

• A simple approximation to the retrospective analysis portion just described is to replace 
the adjoint operator in (19a) by the identity. Since in the current implementation of 
PSAS the forecast error covariance is not dynamically determined, and even with its slowly 
varying forecast error variances it can be thought of as having a time-independent forecast 
error covariance, one might expect that replacing the adjoint by the identity operator in 
(19a) would result in a reasonable retrospective analysis approximation consistent with the 
current underlying PSAS statistics. Todling (2000) has experimented with this idea using 
an identical-twin configuration setup for GEOS and has found a significant improvement 
in the mean error due to lag i — 1 and even to lag i — 2 retrospective analyses. 

(c) The GEOS lag-1 retrospective-based iterated analysis 

When the system is nonlinear, the idea to feed the filter estimate back into the analysis 
equation is particularly attractive, since we expect the filter analysis to be a better estimate of 
the state of the system than the first-guess provided by the model. Indeed, filtering strategies 
making use of such feedback procedures are commonly found in the literature. For instance, 

J azwinski (1970, Theorem 8.2) introduces the so-called iterated EKF, which is suitable for non- 
linear observation operators. Cohn (1997) proposes a similar procedure as an extension to PSAS 
for such operators. Iterative procedures aimed at dealing with nonlinearities of the observation 
operator are sometimes referred to as locally-iterated methods, since the iterations are per- 
formed at a single time. Jazwinski (1970, Theorem 8.3) also presents an iterative procedure 
that is aimed at correcting errors due to the dynamical linearizations required by the EKF. 
This procedure involves integrating the model with a newly estimated trajectory at each iter- 
ation and for this reason it resembles a smoother procedure referred to as the iterated linear 
filter-smoother algorithm. Combining ideas of filtering and smoothing leads to the possibil- 
ity of developing globally-iterated procedures in which the filter analyses may be revised by a 



backward-filter integration within a certain time interval. Most of these iterative procedures are 
inspired by Newton-type methods for solving systems of nonlinear equations (see Navon and 
Legler 1987, and Zou et al. 1993, for reviews of Newton-type methods). 

Motivated by these methods we introduce here a procedure to use the retrospective analysis 
to improve the overall GEOS IAU-based assimilation. At first, the algorithm is based only 
on the lag-1 retrospective analyses. At any given time £&, when a lag-1 retrospective analysis 
w £|^ + i(p) is available we can construct a model-space lag-1 IAU retrospective increment as 

S y a k\k+A a ) = n+ [ w Wi(^)] -y b k\k-i( a ) > ( 23 ) 

which is similar to (14), but is constructed using observations one lag ahead of time £*,. This lag- 
1 retrospective increment can now be used to integrate the GCM over an IAU integration period 
already covered before. This is illustrated schematically in Fig. 2. The diagram resembles the 
regular IAU procedure presented before in Fig. 1. In fact, the top part of the diagram, above 
the horizontal dotted line, is identical to the regular IAU procedure. However, now at, say, time 
t = 12Z we calculate a retrospective analysis by first integrating the transformed PSAS solution 
vector in (18) back in time using the adjoint operation (19a); this is represented in the diagram 
by the southwestward-pointing dashed arrow. A new PSAS-like linear system problem can then 
be solved as in (19b) with the corresponding update (19c), and the lag-1 retrospective analysis 
constructed using (19d), as represented in the diagram by the box tagged “Retro ANA”. In the 
end, a lag-1 retrospective increment at t = 6Z is constructed as indicated in (23), and the GCM 
is integrated for 6 hours using this increment as the tendency term in (9). From this point on, 
the procedure follows the regular IAU schematic until it is time to process the observations at 
t — 18Z when the lag-1 retrospective analysis at t = 12Z can be calculated and the whole cycle 
repeated. The thin blue arrows in Fig. 2 correspond to the retrospective trajectory. In the RIA 
scheme we concentrate on the iterated filter-smoother trajectory represented in the figure by the 
thick solid arrows. At a given analysis time, the relevant iterated PSAS analysis is represented 
in the diagram as the analysis from the lowest PSAS box in a column of the diagram (see thick 
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vertical dashed lines). 


From the diagram in Fig. 2 we see that the retrospective-based iterated analysis results in a 
considerable increase in computational requirements when compared with the regular procedure 
in Fig. 1. ' Each iteration of. the iterated analysis scheme requires one extra 9-hour GCM 
integration and two extra PSAS analyses. Such an increase in the computational cost can only 
be justified if the procedure results in considerably improved analyses. One way to reduce 
the computational burden is by calculating some of the steps in (19) at different resolutions. 
Similarly to the strategy of incremental 4D-var of Courtier et al. (1994), we can for example 
integrate the adjoint of the tangent linear GCM in (19a) at lower resolution than the actual 
model integration (9). Also, the retrospective PSAS-like linear system (19b) can be solved at 
lower resolution than the regular linear system (11) solved in the first step of PSAS. For that 
matter, the calculations in (19a) and (19b) do not even have to be performed at the same 
resolution. This type of approach to reduce computational cost involves the development of 
additional interpolation operators and their corresponding adjoints. 

" Independently of the IAU, in the linear case when the filter portion is actually the Kalman 
filter, it can be shown that to feedback the lag-1 retrospective analysis at, say, tk - i to calculate a 
revised filter analysis at time tk cannot result in an improved filter analysis. In our iterated pro- 
cedure, an optimal analysis could be calculated using the first-guess from the lag-1 retrospective 
analysis if the cross-covariance between the revised first-guess and the observations were prop- 
erly taken into account. In fact, since the retrospective-based iterated analysis procedure here 
amounts to a modified filtering procedure, the optimal gains in this case are similar to the usual 
modified filter gains when the forecast and observations are correlated (e.g., Jazwinski 1970, 
Example 7.5). Since in practice it would be quite difficult to calculate this cross-covariance, we 
choose to neglect the cross-covariance terms all together. 
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4 GEOS experimental results 


fa) Configuration and experimental -setup 

The retrospective analysis procedures of the previous section were implemented as an ex- 
tension to GEOS DAS. The retrospective portion of the GEOS software is compatible with the 
first operational version of GEOS DAS, designed to support NASA’s Earth Observing System 
mission and its Terra satellite. We refer to this earlier operational version as GEOS-3' to avoid 
possible confusion with the considerably upgraded version of GEOS-3 operational at the time of 
this writing. The GEOS-3' GCM operates at a resolution of 1° latitude by 1° longitude and 48 
vertical sigma levels, with a dynamical core essentially like that of Suarez and Takacs (1995). At 
the synoptic hours, PSAS calculates the analysis at a resolution of 2° latitude by 2.5° longitude 
on 20 pressure levels. Details on the implementation of PSAS can be found in da Silva and Guo 
(1996), Guo et al. (1998), and Larson et al. (1998). As we have mentioned in the previous sec- 
tion, GEOS-3' uses the IAU procedure of Bloom et al. (1996) to generate a time-continuous state 
trajectory referred to as the assimilation. For expediency, the experiments performed for the 
present article used both the GCM and PSAS at the coarse horizontal resolution of 4° latitude 
by 5° longitude; the GCM and PSAS vertical resolutions were kept unchanged. We also simpli- 
fiy the experimental configuration by updating the GCM trajectory needed during the adjoint 
integrations only every 6 hours. Except for sea-wind satellite observations, all observation data 
types used in GEOS-3' are included in our experiments. Conventional observations from ships, 
environmental and drifting buoys, surface stations, winds from pilot balloons, aircraft reports, 
and radiosonde stations are used. Cloud track wind retrievals and TOVS geopotential height 
retrievals are used as well. Furthermore, the Wentz (1997) SSM/I-derived total precipitable 
water retrievals are assimilated, though not through PSAS but rather by using the method of 
Hou et al. (2000). 


21 



Four new components are required to implement the retrospective capability in GEOS DAS: 

the adjoint of the tangent linear GCM; the additional PSAS-like operators involved in (19b); 

the linear operator (22a) taking model-space variables into analysis-space variables; and the 

* 

linear operator (22b) taking analysis-space variables into model-space variables. Presently, the 
adjoint of the GCM includes only the hydrodynamics adjoint and the adjoint of a simple diffu- 
sion scheme. Most modifications required to PSAS were quite simple since they only required 
rearranging operators already available in the original PSAS software. Some effort was devoted 
to derive the proper tangent linear and adjoint operators for the transformations (10) and (13), 
because we strove to make sure that the back and forth operations would render minimal error. 
Some of this'work was done by hand, and some was done using the automatic differentiation 
tool of Giering and Kaminski (1998). 

In the present article, only results for the lag-1 (6-hour) retrospective analysis are discussed. 
We compare the results of three experiments conducted over the month of January 1998. To 
minimize possible differences due to spin-up issues, the experiments are actually started on 14 
December 1997, but the results are ignored during this half-month period. Our first experiment 
is taken as the control and it uses the reduced resolution GEOS-3 7 data assimilation system 
mentioned above. The control is referred to as the CTL experiment. In the second experiment, 
referred to as the RA experiment, we also calculate lag-1 (6-hour) retrospective analyses for the 
entire month of January 1998. Since there is no feedback in this experiment, it uses the same 
background fields and OMF time series of the control experiment. The third experiment is aimed 
at evaluating the lag-1 (6-hour) retrospective-based iterated analysis procedure introduced in 
the previous section, and is referred to as the RIA experiment. 

We evaluate the RA and RIA experiments mainly by examining the time-series statistics 
of their corresponding residuals. That is, depending on the case, we calculate time root-mean- 
square (RMS) bias and standard deviation from the differences of the observations with either 
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the forecast, or the analysis, or the retrospective analysis, or the retrospective forecast (see 
below). To ease comparisons, we grid the residuals over 4° latitude by 5° longitude cells on the 
20 pressure levels of the analysis space before calculating any statistics. We calculate statistics 
only from grid-boxes containing 15 or more reports during the month. In the discussion that 
follows, we concentrate on results obtained in the troposphere. 

(b ) Evaluation of the 6-hour retrospective analysis 

We start by comparing the results of the CTL and RA experiments using the set of obser- 
vations assimilated in the CTL experiment. If the 6-hour retrospective analyses are indeed an 
improvement over the regular control analyses we should see that in some mean sense the RA 
observation-minus-analysis (OMA) residuals are reduced in comparison to the OMA residuals 
of the control experiment. As a matter of fact, one can show that in the linear optimal case, 
i.e., when the filter gain is the Kalman gain, 

£{(w° k - H fc w£, fc+1 )(w2 - Hfcw£ |fc+1 ) T } < £{(w£ - Hfcw£| fc ) (w£ - H fc w^ fc ) T } . (24) 

Although there is no guarantee of this holding in general for the suboptimal nonlinear case under 
study, we would like to examine the extent to which it does. In practice, short of perturbing 
the observational data, to assess this quantity we must make the usual ergodic assumption and 
replace the expectation by a time average. Examination of the time RMS biases and standard 
deviations of the OMA residuals when the analyses are either the regular filter analyses of the 
CTL experiment or the lag-1 retrospective analyses of the RA experiment has shown them to 
be virtually identical (results not shown). Therefore, from this point of view we might be led to 
think that there is no payoff in calculating lag-1 retrospective analyses. 

Another way of comparing the quality of two sets of analyses is to compare the skill of 
forecasts issued from them. We expect forecasts issued from retrospective analyses to be superior 
to regular forecasts for at least their total lag period, 6 hours in the lag-1 case here, since their 
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initial conditions have had the benefit of observations that far into each forecast. Hence, we 
compare the OMF residual statistics of the regular filter forecasts of the CTL experiment and of 
the so-called retrospective forecasts issued from the lag-1 retrospective analyses. Since the OMF 
residuals from a regular GEOS DAS run, such as the CTL experiment, involve 6-hour forecasts 
that are produced from partly integrating the GCM with the IAU forcing for 3 hours and partly 
integrating the GCM for another 3 hours without the influence of the IAU tendencies (see Fig. 
1), we must use the retrospective analyses carefully when constructing OMF residuals from them. 
To make a fair comparison, we calculate OMF residuals from the 6-hour retrospective analyses 
following a forecasting procedure based on IAU. For each available retrospective analysis for the 
entire month of January 1998 a retrospective forecast is issued following the schematic shown in 
Fig. 3. As illustrated in the figure, the retrospective OMF residuals at, say, 12Z are calculated by 
converting the 6Z retrospective analysis to the model space and constructing the corresponding 
increment on the model space, following (23). This retrospective analysis increment is used as 
a tendency term during a 6-hour GCM integration, started at 3Z. At the end of the 6-hour 
integration the retrospective tendency term is turned off, by setting a = 0 in (9), and the model 
is left to run free for another 3 hours, after which the OMF residuals at 12Z can be calculated 
using the observations at that time. 

Using these retrospective forecasts, Fig. 4 shows the time RMS bias (top panels) and stan- 
dard deviation (bottom panels) for the radiosonde geopotential height OMF residuals for the 
CTL (solid curves) and RA (dashed curves) experiments averaged over the western (left) and 
eastern (right) quadrants of the Northern Hemisphere, for latitudes higher than 20N. These 
two domains are chosen because they represent the largest concentration of radiosondes over 
the globe. We see from the top panels that, in the RMS bias sense, the forecasts from the 
lag-1 retrospective analyses correspond to a considerable improvement over the regular GEOS 
DAS analyses. However, the bottom-left panel shows that for the RMS standard deviations the 
retrospective forecasts are considerably degraded compared to the regular forecasts over what 
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is mostly North America; results are roughly neutral over most of Europe and Asia, as seen 
from the bottom-right panel. Figure 5 displays similar quantities but now for the zonal wind 
radiosonde OMF residuals. Except in the Northwestern region between pressure levels 700 hPa 
and 400 hPa, the RMS bias for the zonal wind radiosonde OMF residuals is improved when fore- 
casts are issued from the lag-1 retrospective analyses. In this same region, the zonal wind OMF 
standard deviation [panel (b.l)] shows a minor deterioration at levels below 400 hPa, much less 
than that seen in the OMF heights in Fig. 4b. 1; minor improvement in the standard deviations 
are seen above 400 hPa. Over the Northeastern region a minor but consistent improvement is 
observed in both the RMS bias and standard deviation, as indicated in the panels on the right. 

The statistics of OMF residuals for other variables and other observing systems can also be 
examined. Figure 6 shows the time RMS biases (top panels) and standard deviations (bottom 
panels) for the TOVS geopotential height OMF residuals. Since TOYS provides global coverage 
in the course of a single day, the spatial averages now cover the entire Northern Hemisphere 
(left panels) and Southern Hemisphere (right panels). We see considerable improvement in the 
OMF biases and standard deviations from the retrospective forecast residuals. Interestingly, 
the standard deviation results over the Northern Hemisphere [panel (b.l)] contradict the dete- 
rioration observed in the radiosonde geopotential height OMF residuals [panel (b.l) of Fig. 4]; 
We attribute this contradiction over North America to contradictions between the geopotential 
height observations from the radiosondes and the TOVS retrievals themselves and not to the 
retrospective analysis procedure. 

Another quantity we have studied is simply the spatial average of the residuals time mean. 
Though we expect considerable cancellation of errors in this quantity, it still serves as an indicator 
of the overall behavior of the residuals and of the underlying procedure used to produce them. 
Figure 7 shows the time mean OMF residuals for the CTL 6-hour first-guesses (forecasts) and 
the lag-1 retrospective forecasts. The globally-averaged time means for TOVS and radiosonde 
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geopotential height OMF residuals are displayed in panels (a) and (b), respectively. We see 
mostly a reduction in the time mean residuals when the retrospective forecasts are used instead 
of the regular forecasts, with some overshooting at levels below 700 hPa for the TOVS residuals. 
The zonal and meridional wind components of the radiosonde OMF residuals are displayed in 
panels (c) and (d), respectively, and again we see an overall reduction when the retrospective 
forecasts are used, with some overshooting of the mean meridional wind around 150 hPa. 

In terms of the metrics presented here for the nonlinear suboptimal case of the GEOS applica- 
tion, we see clear benefit in producing the 6-hour forecasts from the lag-1 retrospective analyses 
over the regular GEOS DAS forecasts. This serves to indicate improved analysis quality with 
the RA scheme. This also serves as further motivation to consider the iterated retrospective 
analysis procedure proposed in the previous section, since it makes direct use of the retrospective 
forecasts (see Fig. 2). 

(c) Evaluation of the 6-hour retrospective-based iterated analysis 

We now evaluate the performance of the 6-hour (lag-1) retrospective-based iterated analysis 
scheme. We start by comparing the OMA residuals between the CTL and the RIA experiments. 
Figure 8 shows the globally- averaged time RMS bias for TOVS and radiosonde geopotential 
height OMA residuals [panels (a) and (b), respectively], and for the zonal and meridional com- 
ponents of the radiosonde winds [panels (c) and (d), respectively]. Although small, we now 
actually see improvement in the OMA residuals due to the iterated analysis. To the extent that 
the expectation can be replaced by the RMS time mean, the inequality (24) holds when w^ +1 
corresponds to the iterated analysis, at least in a globally-averaged sense. We have examined 
the time RMS standard deviations of the OMA residulas of both the CTL and RIA experiments 
and have found it to change very insignificantly. 

Even though small, the improvement due to the RIA scheme seen in Fig. 8 is also visible 
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directly from the time series of the globally-averaged OMA residual RMS bias. Furthermore, 
this improvement is seen not only for TOVS and radiosonde OMA residuals, but for other 
instruments as well. An illustration is presented in Figure 9 by displaying the globally-averaged 
RMS bias of the zonal (top) and meridional (bottom) cloud-track wind OMA residuals at 200 
hPa. The thin curves correspond to the OMA residual time series from the CTL experiment 
and the thick curves are for the RIA experiment. The global reduction in the RMS bias can be 
as much as 1 m s -1 at times. This confirms the reduction in the globally-averaged time RMS 
bias of the radiosonde OMA residuals observed in panels (c) and (d) of Fig. 8 around the same 
pressure level. 

Frequently, changes made to assimilation systems are evaluated and validated by making 
comparisons with independent observations, that is, observations which are not assimilated 
by the system. Data withholding experiments are commonly used to assess the impact of a 
particular observing system and can also be used to evaluate the impact of system changes (e.g., 
Bouttier and Kelly 2001, and references therein). Here we choose to validate the change in 
the 200 hPa winds of Figs. 8 and 9 by using wind observations from the Global Aircraft Data 
Set (GADS) of the British Airways Boeing 747-400 flights, and by using further aircraft wind 
observations from the Aircraft Communications, Addressing, and Reporting System (ACARS). 
Neither of these observation types were used in our assimilation experiments and therefore they 
provide independent checks. 

The GADS wind observations have been shown by Rukhovets et al. (1998) to be of value 
to GEOS DAS if used regularly in the PSAS analyses. This suggests that any changes made to 
GEOS that show its analyses to draw more closely to these observations, even when they are 
not assimilated, should be considered an improvement. With that in mind, we used a dataset 
of the January 1998 GADS observations to construct OMA residuals for the analyses of both 
the CTL and RIA experiments. Figure 10 shows maps of the time RMS standard deviation 
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of the zonal (top) and meridional (bottom) GADS winds OMA gridded residuals for the CTL 
experiment subtracted from the same quantity for the RIA experiment (RIA-minus-CTL). The 
color scheme in the figure indicates that blue (negative values) corresponds to improvements 
due to the RIA procedure. Though we see areas where the impact of RIA is neutral or negative, 
in most places the GADS observations are closer to the analyses of the RIA experiment than to 
those of the CTL experiment. 

Similarly, Fig. 11 shows the differences of Fig. 10, but now for the ACARS wind OMA 
residuals. The maps are focused over North America since that is where the majority of the 
observations are concentrated in this case. Relatively neutral results are seen in the meridional 
component of the wind (bottom map), but undeniable improvement due to the RIA scheme is 
seen in the zonal component of the wind (top map). 

Ultimately, as emphasized by CST94, one of the main motivations for performing retrospec- 
tive analysis is to produce the best possible dataset for climate research. As such, it is important 
to examine the climatological impact of changes induced by the RIA procedure. Since the results 
of the experiments discussed here are still preliminary we do not want to dwell too much on the 
significance of performing RA and RIA for the purposes of improving the climatological aspects 
of the assimilation strategy — recall that our experiments are for a very low resolution version 
of GEOS-3'. Still, we cannot avoid looking more closely to see what is the climatological impact 
of changes such as those observed in the wind field. In fact, the significance of the RIA impact 
on the upper-level winds can be seen more clearly by looking directly at the monthly-averaged 
winds. For instance, Fig. 12 shows the zonally- averaged January 1998 monthly mean merid- 
ional wind (top) for the RIA experiment and its difference from the CTL experiment (bottom). 
The bottom panel shows a distinct tropical wind strengthening at the upper levels and a slight 
weakening at the mid- to lower levels when the RIA scheme is used. 

This change in the tropical meridional wind affects the Hadley circulation. To see the 
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meridional circulation, we calculate the mass stream function ip by integrating the zonally- 
averaged monthly mean meridional wind using the expression 


iP = 


2i tR cos < 


H dp ' , 


(25) 


' Ptop 


where v stands for the meridional wind, the operator • represents the time mean, the operator 
[•] represents for the zonal average, R is the mean radius of the earth, g is the gravity constant, <j> 
is the latitude, and the integral is from p top = 10 hPa to a pressure p down to the surface. Figure 
13 shows the January 1998 mass stream function for both the CTL (top) and RIA (bottom) 
experiments. We see a clear enhancement of the Hadley circulation when the RIA procedure is 
used, with the mass stream function peaking at about 16 X 10 10 kg s~ 1 in contrast to the weaker 
peak of 12 x 10 10 kg s -1 for the circulation of the CTL experiment. Although we do not expect 
the circulation to be well-represented at the coarse resolution we use in our experiments here, it 
is much closer to the circulation pattern of the full-resolution, 1° latitude by 1° longitude GEOS 
DAS (not shown), with its tropical circulation peaking at 18 X lO 10 kg s -1 . This suggests that 
the RIA scheme has the potential for improving climatologically relevant features. 


Finally, we compare the skill of 5-day forecasts issued from the CTL and the RIA analyzed 
fields. These are initialized as in Fig. 3. Since our experiments are confined to the single month of 
January 1998, we have few independent samples for this comparison. We issued 5-day forecasts 
starting from 2 January 1998 every 3 days until 26 January 1998, to have a small sample of 9 
5-day forecasts. We verified that the overall conclusions and skills calculated from this small 
ensemble were not affected by the size of the sample by reducing the size of the sample to 5 
members and performing cross-validation. As a measure of forecast skill we calculated anomaly 
correlations and RMS errors [e.g., von Storch and Zwiers 1999, Eqs. (18.17) and (18.18)]. Both 
the CTL and RIA forecasts were verified against their own analyses. Anomalies were calculated 
using a 10-year climatology obtained from ECMWF operational analyses for the period of 1988 
to 1997, and interpolated to the resolution of our experiments. We should say that the scores 
shown below are not representative of the scores of the operational GEOS-3 data assimilation 
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system, which operates at higher resolution than that of the experiments here. 

Figure 14 shows the anomaly correlations for the 500 hPa geopotential height field calculated 
over four different regions for the 5-day forecasts issued from the CTL (solid curves) and RIA 
(dashed curves) analyses. We see that over the Northern Hemisphere extratropics (top-left 
panel) forecasts from RIA analyses are of similar skill as forecasts from the control analyses, 
at least up to day 4. Over North America (bottom-left panel) the forecast skill from RIA 
analyses show some deterioration when compared against the skill of the CTL forecasts. As when 
studying the OMF residuals obtained from the retrospective forecasts using the RA analyses, 
this deterioration over North America might be related to contradictions in the observing system 
over this area (see Figs. 4b. 1, 6b. 1). In fact, this seems to be an issue confined to this region 
since, for example, over the Southern Hemisphere (top-right panel) and Europe (bottom-right 
panel) we see improvement in skill when the 5-day forecasts are issued from the RIA analyses. 

As a final illustration comparing the 5-day forecast skill from the CTL and RIA experiments 
we examine the RMS error of the tropical wind fields at 850 hPa and 200 hPa. Figure 15 displays 
these quantities for both the zonal (left panels) and meridional (right panels) components of the 
wind. The RMS errors at 850 hPa are virtually identical, while at 200 hPa we see a slight 
improvement when using forecasts from the analyses of the RIA experiment. Although these 
are small improvements they serve as further confirmation of what we have seen previously when 
comparing the analyses of the CTL and the RIA experiments with independent observations. 

5 Conclusions 

A central purpose of atmospheric data assimilation is to produce the best possible estimate 
of the the state of the atmosphere at any single time. In theory this can be accomplished by 
using smoothing techniques since they are aimed at maximizing data usage through inclusion of 
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observations in the past, present, and future of the time an estimate is sought. In the context 
of sequential data assimilation, the fixed-lag Kalman smoother (FLKS) provides a particularly 
attractive framework. The FLKS formulation is based fully on the underlying filtering strategy. 
Its standard formulation requires no error covariance information beyond what is required by the 
filtering approach. Indeed, the FLKS can be separated into a filter portion and a retrospective 
analysis (RA) portion and this separation renders practical implementation of FLKS-based 
procedures a relatively simple extension of an already existing (filter) analysis scheme. 

Two different types of retrospective procedures are investigated in the present work. The 
first is the original FLKS-based formulation referred to simply as RA. The second is an iterated 
version of the original algorithm, referred to as RIA, in which lag-1 retrospective analyses are 
used to revise the previously calculated filter analysis. Both these procedures are implemented 
as an extension of the Goddard Earth Observing System (GEOS) Data Assimilation System 
(DAS). The new components required for implementing a retrospective capability in GEOS 
DAS are the adjoint of the tangent linear model of the GEOS general circulation model (GCM); 
the rearrangement of a few operators already available in the physical-space statistical analysis 
system (PSAS) of GEOS DAS; and the development of the tangent linear and adjoint operators 
responsible for transforming between model-space variables and analysis-space variables as well 
as the adjoint of the operator transforming analysis-space variables into observables. The adjoint 
of the tangent linear GCM used in the present work includes the adjoint of the tangent linear 
hydrodynamics and the adjoint of a simple diffusion term; the adjoint of the physics is not 
included. 

Only results for the 6-hour (lag-1) retrospective analysis were studied here. Although close 
examination of the observation-minus- analysis (OMA) residuals seemed to suggest a rather 
neutral benefit from the lag-1 retrospective analysis, we saw improvement in the 6-hour forecasts 
issued from these lag-1 retrospective analyses: these so-called retrospective forecasts are a closer 
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match to the observations than the regular GEOS DAS forecasts. This improved 6-hour forecast 
skill motivated the investigation of the RIA scheme since this scheme makes explicit use of the 
retrospective forecasts. Evaluation of the analyses from the RIA procedure indicated them to be 
closer to the observations than the usual PSAS analyses. The OMA residuals for independent 
observations not used during the assimilation further confirmed some of the improvements due 
to the RIA scheme. More significant improvements were seen when examining climatologically 
important fields such as the mass stream function describing the meridional wind circulation. 
Lastly, anomaly correlations and root-mean-square errors from a small sample of 5-day forecasts 
indicated a mild improvement in skill scores when analyses from the RIA procedure were used 
for the 5-day forecasts instead of the regular GEOS analyses. Although the skill scores were not 
improved everywhere over the globe, they were improved generally. 

Much work remains to be done to show that retrospective analysis is a worthwhile extension 
to the usual PSAS analysis of GEOS. The present work used only a reduced-resolution version 
of GEOS DAS and a study with a higher resolution version is necessary. Just as in four- 
dimensional variational procedures, one of the main features of the retrospective analysis is 
its capability to incorporate backward-propagated information into the retrospective analyses 
via the model adjoint. This should, in particular, result in improved representation of synoptic 
features. Synoptic evaluation has not been explored in the present work and awaits full-resolution 
experimentation. Furthermore, it is evident that instead of windowing the observations in 6- 
hour batches, as is commonly done in 3D-var systems, much can be gained by assimilating 
observations at their proper time, as is done in 4D-var. An alternative rapid update cycle 
strategy is being considered for this purpose at the DAO and it will be important to investigate 
the impact of retrospective analysis in this context. These studies will be done primarily in the 
context of the newly developed finite- volume data assimilation system. In any case, the results 
presented here show promise for the development of retrospective analysis capabilities in this 
new assimilation system. 
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Appendix A 

Retrospective gains as a function of filter variables only 


The purpose of this appendix is to derive the alternative expression (16) for the retrospective 
gain matrix (see also Zhu et al. 1999). Using (6) and (7) with £ = 1, 2, ..., j we have 


Kfc-ijfc — 


Kfc_2|jfc — 




pa a T xrTin-1 

P{_ 1]k _ 2 [(I - E*_ x ] Hjr^ 1 , (A. la) 

(Pfc-i lfc - 2 |fc- 1 ) T Ar,*- 1 Hjrj 1 

Pfc_2|fc-2-^-i-i,fc-2 [(i ~ Kfc_i|fc_iH^_i) A fcjfc _ 1 ] H fc r fc 

Pfc-2|jk-3 [O' _ Kyt_2|A;-2Hfc-2) T A^_ 1)fc _2] [(I — A^j.^] 

xH^ 1 , ' ( A - lb ) 

H-Xk-i-i [(I-K H |H H H) r AL i+1 ^] ••• [(I-K^n^H^ifAj^] 
xHfr^ 1 . ( A - lc ) 


This can be written generally as in (16) or, making explicit use of (1c) for the filter gain matrix, 
we can also write 


K 


k-t\k 


— r k-l\k-t-l 


II (I - JAC_, 


T-r-1 


uiT 


k 1 - k ’ 


(A.2) 


j—k — £-\- 1 


which shows that, as pointed out in the main text, the retrospective gains depend only on filter 
quantities. Todling et al. (1998) have pointed out that the retrospective portion of the FLKS 
implicitly accounts for model error. The equation above serves to re-emphasize that remark as 
it shows that the retrospective gains depend directly on the forecast error covariance matrix, 
which is the filter quantity containing model error covariance information. 
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Figure 1: Schematic representation of the incremental analysis update (IAU) procedure. The 
dashed arrows represent the 3-hour GCM integration that provides the first-guess (forecast) to 
PSAS. At each analysis time PSAS uses observation-minus-forecast (OMF) residuals to calculate 
an updated state estimate (analysis; vertical dotted lines). The analysis-minus-forecast difference 
is converted to a model-space tendency term used to force the GCM during a 6-hour integration 
around the analysis time; this is the IAU period represented by the solid thick arrows. The 
cycle is repeated with a 3-hour GCM integration, without the IAU tendency term, to provide 
the first-guess for the next analysis time. The line formed by the solid arrows represents a 
time-continuous IAU trajectory, referred to as the assimilation. (Similar to Fig. 1 of Bloom et 
al. 1996). 
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Figure 2: Schematic representation of the lag-1 iterated retrospective data assimilation pro- 
cedure. Dashed north-eastward pointing arrows represent GCM first-guess integration; solid 
eastward-pointing arrows represent GCM integration forced by IAU (thick) and retrospective 
IAU increment (thin). Dashed south-westward-pointing arrows represent 6-hour adjoint model 
integrations. The boxes labeled “Retro ANA” stand for the PSAS application in (19b). The 
retrospective assimilation is used to provide a revised first-guess that is further used to revise 
the filter analysis at each synoptic time. 



Figure 3: Schematic of the procedure to issue forecasts from retrospective analyses using the 
IAU framework. Arrows are similar to those in Fig. 1. The main purpose of the retrospective 
forecast is the calculation of the OMF residuals indicated by the “Retro OMF” box. 
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Figure 4: Time root-mean-square (RMS) bias (panels a) and standard deviation (panels b) for 
the radiosonde geopotential height OMF residuals for the control experiment (solid curves) and 
for forecasts from the lag- 1 retrospective analyses from the RA experiment (dashed curves). 
Panels 1 on the left are for the Northwestern quadrant of the globe defined between longitudes 
180W-0 and between latitudes 20N-90N; panels 2 on the right are for the Northeastern quadrant 
of the globe between longitudes 0-180E and latitudes 20N-90N. Units are in meters on the 
abscissa and hPa on the ordinate. 
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Figure 7: Globally- averaged time mean OMF residuals for the CTL experiment (solid curves) 
and for the retrospective forecasts (dashed curves). Panel (a) is for the geopotential height 
TOVS retrievals residuals; panels (b)-(d) are for the geopotential height, zonal wind, meridional 
wind radiosondes residuals, respectively. 
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Figure 8: Globally- aver aged time RMS bias of observation-minus- analysis (OMA) residuals for 
the CTL experiment (solid curves) and retrospective assimilation from the RIA experiment. 
Panels are arranged as in Fig. 7: panel (a) is for the TOVS geopotential height residuals; panels 
(b)-(d) are for the radiosonde geopotential height, zonal wind, and meridional wind residuals, 
respectively. 
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for the CTL experiment, and the thick curves are for the RIA experiment. Units are m s' 
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Figure 10: Maps of the 200 hPa time RMS standard deviation of the GADS winds OMA residuals 
of the CTL experiment subtracted from the same quantity for the RIA experiment. The top 
map is for the zonal wind and the bottom map is for the meridional wind. Units are in m s -1 . 
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Figure 11: Same as in Fig. 10, but for the ACARS wind OMA residuals. Only North America 
is displayed since it corresponds to the area where the bulk of these observations are. 





Figure 13: Mass stream function for CTL (top) and RIA (bottom) experiments. Units are i 


Figure 14: Anomaly correlations for 500 hPa geopotential heights for CTL (solid curves) and 
RIA (dashed curves) experiments. 














